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ABSTRACT 


Q-ball  imaging  (QBI)  is  a  high  angular  resolution  diffusion  imaging  (HARDI)  technique  which  has  been 
proven  very  successful  in  resolving  multiple  intravoxel  fiber  orientations  in  MR  images.  The  standard 
computation  of  the  orientation  distribution  function  (ODF,  the  probability  of  diffusion  in  a  given 
direction)  from  q-ball  data  uses  linear  radial  projection,  neglecting  the  change  in  the  volume  element 
along  each  direction.  This  results  in  spherical  distributions  that  are  different  from  the  true  ODFs.  For 
instance,  they  are  neither  normalized  nor  as  sharp  as  expected,  and  generally  require  post-processing, 
such  as  artificial  sharpening  or  spherical  deconvolution.  In  this  paper,  a  new  technique  is  proposed  that, 
by  considering  the  solid  angle  factor,  uses  the  mathematically  correct  definition  of  the  ODF  and  results  in 
a  dimensionless  and  normalized  ODF  expression.  Our  model  is  flexible  enough  so  ODFs  can  either  be 
estimated  from  single  q-shell  datasets,  or  exploit  the  greater  information  available  from  multiple  q-shell 
acquisitions.  We  show  that  this  can  be  achieved  by  using  a  more  accurate  multi-exponential  model  for  the 
diffusion  signal.  The  improved  performance  of  the  proposed  method  is  demonstrated  on  artificial  data  and 
real  HARDI  volumes. 


Key  words:  Orientation  distribution  function  (ODF),  q-ball  imaging  (QBI),  high  angular  resolution 

diffusion  imaging  (HARDI),  solid  angle. 
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INTRODUCTION 


Diffusion-weighted  magnetic  resonance  imaging  (DWMRI)  provides  valuable  information  about  the  fiber 
architecture  of  neural  tissue  by  measuring  the  diffusion  of  water  molecules  in  three-dimensional  (3D) 
space.  The  microscopic  diffusion  function  may  be  measured  by  using  the  model-free  diffusion  spectrum 
imaging  (DSI)  (1),  which  uses  the  direct  Fourier  inversion  of  the  diffusion  signal.  This  technique  is  time 
intensive,  as  it  measures  the  diffusion  signal  on  a  3D  (e.g.,  11x11x11)  Cartesian  lattice.  Thus,  an 
alternative  approach  based  on  sampling  only  on  one  or  multiple  spherical  shells  in  q-space  has  been 
proposed,  referred  to  as  high  angular  resolution  diffusion  imaging  (HARDI)  (2).  The  spherical  shell, 
being  a  2D  manifold,  includes  a  number  of  measurement  points  which  grows  quadratically  with  the 
desired  angular  resolution,  as  opposed  to  cubically  with  the  spatial  resolution  in  the  entire  3D  lattice  of 
DSI. 

While  the  3D  probability  density  function  (PDF)  of  the  diffusion  is  helpful  in  studying  the  tissue 
microstructure,  the  orientation  distribution  function  (ODF)  -  the  marginal  probability  of  diffusion  in  a 
given  direction  -  is  the  quantity  of  interest  for  mapping  the  orientation  architecture  of  the  tissue.  Q-ball 
imaging  (QBI),  (3),  is  a  widely  used  acquisition  scheme  for  HARDI,  from  which  ODFs  can  be 
reconstructed  through  a  spherical  tomographic  inversion  called  the  Funk-Radon  transform.  This 
technique’s  simplicity  and  its  ability  to  resolve  intravoxel  fiber  orientations  have  made  it  popular  for  fiber 
tracking  and  characterizing  white  matter  architecture.  A  number  of  recently  proposed  methods  have 
turned  QBI  into  a  very  efficient  and  robust  technique  (4)-(6).  Moreover,  a  few  works  have  suggested 
exploiting  data  from  multiple  q-shells  to  benefit  from  the  high  signal-to-noise  ratio  (SNR)  and  high 
angular  contrast-to-noise  ratio  (CNR)  of  the  data  acquired  at  respectively  low  and  high  b-values,  (3),  (7)- 
(9).  Using  multiple  q-shells  also  allows  us  to  employ  richer  models  for  the  diffusion  signal,  as  discussed 
in  this  paper. 

However,  the  definition  of  the  ODF  used  in  the  original  QBI  is  different  from  the  actual  marginal  PDF 
of  diffusion  in  a  constant  solid  angle.  It  is  computed  as  a  linear  radial  projection  of  the  PDF,  which  does 
not  take  into  account  the  quadratic  growth  of  the  volume  element  with  respect  to  its  distance  from  the 
origin  (see  the  “General  ODF  Definition”  section  and  Fig.  1  for  more  details).  This  inaccurate  formulation 
generally  distorts  the  ODF,  produces  non-distribution  functions,  and  has  created  the  need  for  artificial 
post-processing  such  as  manual  normalization  and  sharpening. 

In  this  paper,  we  re-derive  the  ODF  expression  for  QBI  via  Fourier  analysis,  this  time  starting  from  the 
proper  definition  of  the  ODF  in  constant  solid  angle.  We  show  that  this  results  in  an  inherently 
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normalized  and  dimensionless  expression.  In  addition,  we  illustrate  through  our  experiments  that  the 
ODFs  are  naturally  sharp  and  that  multiple  fiber  orientations  are  thus  better  resolved.  We  also  provide  a 
general  formulation  for  multiple  q-shell  QBI,  and  demonstrate  the  improvement  achieved  by  considering 
the  information  from  multiple  q-shells  and  using  richer  multi-exponential  models.  Furthermore,  by 
making  use  of  the  spherical  harmonic  basis,  we  demonstrate  that  the  implementation  of  the  new, 
mathematically  correct  expression  is  as  straightforward  as  that  of  the  original  formula,  or  maybe  even 
simpler,  considering  that  further  sharpening  (post-processing)  is  not  necessary. 

This  paper  extends  our  previous  conference  versions  for  single  (10)  and  multiple  q-shells  (11).  In 
particular,  we  provide  complete  mathematical  proofs,  a  regularization  scheme,  and  additional  validation 
and  comparisons.^ 


METHODS 

General  ODF  Definition 

The  PDF  of  the  diffusion  of  water  molecules,  P{f),  gives  the  displacement  probability  P{f)dv  of  a 
molecule,  initially  placed  at  the  origin,  to  be  in  the  infinitesimal  volume  dv  located  at  f  after  a  certain 
amount  of  time.  We  assume  this  function  to  be  symmetric  (i.e.  p{-  r)  -  P{r)),  which  is  a  quite  common 

assumption  in  DWMRI.  The  PDF  can  be  presented  in  Cartesian  coordinates  with  f  =  {x,  y,  zf  and 
dv  =  dxdydz  ■  However,  for  mapping  the  orientation  architecture  of  the  tissue,  the  representation  which 
mostly  interests  us  is  in  the  standard  spherical  coordinates,  parameterized  by  where  f  -ru,  with 

u{e,(t>)  =  (sin  ^  cos  sin  ^  sin  cos  (9)^  being  the  unit  direction  vector.  The  volume  element  in  this  case  is 
dv  =  r^drdO.  with  dQ.  =  ?,m6d6d^  being  the  infinitesimal  solid  angle  element. 

^  After  our  conference  paper  was  accepted  and  its  extension  to  multiple  shells  was  submitted,  a  parallel 
and  independent  work  was  published  (14),  where  the  proper  definition  of  the  ODF  was  considered  in 
single  q-shell  QBI.  However,  in  addition  to  not  considering  multiple  shells  and  the  richer  model  as  done 
here,  the  authors  of  (14)  take  the  integral  of  the  diffusion  signal  on  a  circle  and  not  on  the  entire  plane, 
and  that  results  in  a  different  formula  which  is  not  necessarily  normalized  and  leads  to  other  potential 
inaccuracies.  (See  the  “Q-ball  Imaging  ODF  Reconstruction”  section  for  further  details  and  comparison). 
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We  denote  by  ODF{u)dQ.  the  probability  of  diffusion  in  the  direction  u  through  the  solid  angle  d^, 
which  can  be  computed  by  integrating  the  displacement  probabilities,  i.e.,  P{f)dv  =  P{ru)r^drdQ. ,  for  all 
magnitude  r ,  while  keeping  u  constant: 

ODF{u)dQ.=  \  P(ru)r^ drdQ. , 

Jr=0 

or  simply: 

ODF{u)  =  ^°°^P{ru)r^dr.  [1] 

The  above  definition,  which  is  normalized  and  dimensionless,  is  the  integral  of  the  probability  values 
in  a  cone  of  “very  small”  constant  solid  angle  (Fig.  1(a)).  This  correct  definition  was  used  for  instance  by 
the  authors  of  (1)  in  DSI,  where  P{r)  was  first  computed  from  the  diffusion  data  via  Fourier  inversion 
and  then  integrated  to  calculate  the  ODF,  and  also  in  (12)-(13)  for  diffusion  tensor  imaging  (DTI),  where 
the  ODF  was  analytically  computed.  However,  the  original  expression  for  ODF  reconstruction  in  HARDI, 
and  specifically  QBI  (3),  is  different  from  Eq.  [1],  in  the  sense  that  the  integral  is  not  weighted  by  the 
important  (and  mathematically  correct)  factor  (Fig.  1(b)).  To  the  best  of  our  knowledge,  the  only  paper 
which  has  so  far  considered  this  factor  in  (single  shell)  QBI,  is  a  very  recent  parallel  work  (14)  (published 
independently  after  a  conference  version  of  our  paper  (10)  had  just  been  accepted),  where  the  ODF  is 
approximated  from  the  q-sheh  using  Eq.  [1].  (See  the  “Q-bah  Imaging  ODE  Reconstruction”  section  for 
details  and  comparison.) 

Computing  the  ODE  without  the  factor  would  be  equivalent  to  assuming  the  PDE  to  be  F(r)/|fp  ,  as 

P{ru)dr  =  .  This  radial  projection  gives  an  artificial  weight  to  plf)  which  is,  respectively, 

Jo  ^  Jo 

too  large  and  too  small  for  points  close  to  and  far  from  the  origin,  and  in  fact,  the  computed  quantity 
would  be  different  just  as  the  zeroth  moment  of  a  one-dimensional  function  p{r):-  P{ru)  is  different  from 

its  second  moment.  Eor  instance,  a  consequence  of  not  including  the  factor  is  that  the  computed  ODE 
will  not  be  necessarily  normalized,  and  an  artificial  normalization  factor  will  be  required.  Moreover,  the 
ODE  will  not  be  dimensionless,  since,  given  that  p(f)  has  the  dimension  of  Lr’  (L  being  the  length 

dimension),  the  dimensions  of  p{f)r^dr  and  P{f)dr  are  respectively  1  and  L“^. 
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As  an  example  intended  for  comparison,  we  compute  the  ODFs  with  and  without  in  the  case  of  DTI, 
with  the  following  standard  Gaussian  PDF: 


where  D  is  the  covariance  matrix  (proportional  to  the  diffusion  tensor).  The  computed  ODFs  are: 

ODFiu)^ - ^ 

with  \  \—  /  T  ^  \~ 

A7r\D\^[u^D-^uf 
ODF  {u)=- - - r’ 

without  T,  I  if  /  T  ,  \  — 

\K\DY\ffD-^uY 

where  Z  is  the  normalization  constant  that  subsequently  needs  to  he  computed  and  considered  in 
ODF  {u)  (see  (3)).  An  example  of  this  pair  of  ODFs  is  illustrated  in  Fig.  2.  (No  min-max  normalization  is 

without 

used  in  any  of  the  figures.) 

Next,  we  derive  a  closed-form  expression  for  the  ODF  in  QBI  using  the  correct  -weighted  integral. 

Q-ball  Imaging  ODF  Reconstruction 

Let  E{q)  he  the  3D  Fourier  transform  of  P{f).  We  have  the  values  of  E{q)  measured  on  a  q-hall,  i.e.,  the 
frequencies  with  constant  norm  |^|  =  ^g,  as  e{u):=  E{q(,u)  =  ,  where  S{u)  is  the  HARDI  signal  and 

Sq  is  the  non  diffusion-weighted  (or  BO)  image.  In  addition,  since  the  diffusion  signal  at  ^  =  0  is  5^ ,  one 
can  see  that  £(0)  =  !.  Alternatively,  £(0)  is  the  zero  frequency  of  a  PDF  which  is  its  integral  over  the 
entire  space,  yielding  1 . 

Our  mathematical  derivation  is  based  on  the  following  two  fundamental  facts  from  Fourier  analysis: 

•  The  Fourier  transform  of  p(r)|Fp  is  -V^E{q),  where  is  the  Laplacian  operator  (proof 
presented  in  Appendix  A). 
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•  For  a  symmetric  function  with  the  3D  Fourier  transform  function  f{q),  and  for  the 

arbitrary  unit  vector  u,  we  have  that  V"  f{^fu)dr  =  -^\\  f{q)d^q,  where  is  the  plane 

Jo  JJ.± 

perpendicular  to  u  (proof  presented  in  Appendix  B). 

Combining  these  statements  with  Eq.  [1]  leads  to 


ODF{u)  =  -^\lXE{q)d^q 


Now,  without  loss  of  generality,  we  choose  our  coordinates  such  that  z-u,  thus  making  u''~  the  qx-qy 
plane.  We  then  use  the  following  expansion  for  the  Laplacian  in  spherical  coordinates  {q,0,^): 


V^E{q)^~iqE)  +  \vlE, 


q  dq 


where  is  the  Laplace-Beltrami  operator  which  is  defined  independently  of  the  radial  component  q,  as 


=  — — 

*  sin^d^ 


0  =  —  and  using  the  expression  d^q  =  qdqdc/)  for  the  surface  element,  which  yields 
2 


sin^ 


dE 


1  d^E 


do )  sin  0  d(p 


.  The  surface  integral  on  the  qx-qy  plane  is  computed  by  fixing 


I  C  '2.71  /I 

ODF(z)  =  -^£  £  V  E{q)qdqd(p 


1 

1  /•2;r  1  8^  ^  X  1  9  ^ 


qdqdcp 


We  can  see  that  the  integral  of  the  first  term  is  constant  and  independent  of  E{q)  and  its  derivatives: 


dq^ 

-[E  +  qE^l 
Jo 

^E[o.)-E[0)+[dEX__A‘iEX., 

=  -l 


qdqd(l>  =  —271  ■ 


7 


where  the  subscript  indicates  the  partial  derivative.  We  assumed  (the  standard  assumptions)  that  the 
diffusion  signal  and  its  radial  derivative  go  to  zero  (sufficiently  fast)  as  q  —^oo ,  and  also  that  the 
derivative  is  hounded  at  the  origin.  Therefore  we  have 

ODF{z)  =  - - ^-^C"r-VlE{q)dqd(p,  [2] 

^  ^  4;r  8;r"-'o  q  ”  ^  ^ 

TT 

while  ^  =  —  is  kept  constant  in  the  integration. 

2 

To  compute  the  integral  of  the  second  term,  the  values  of  E{q)  are  required  in  the  entire  q-space.  The 
above  equation  could  be  used  for  example  in  the  DSI  modality,  where  direct  computation  of  the  ODF 
from  the  diffusion  signal  would  eliminate  the  need  for  3D  Fourier  inversion.  In  QBI,  however,  the  values 
of  E{q)  are  not  available  in  the  entire  q-space.  Thus,  we  need  to  approximate  E{q)  from  the  values 
measured  on  the  q-ball.  In  this  work,  we  consider  the  following  radial  mono-exponential  model: 

^2  ^2 

E{q)  =  E{qji)^  =  E{u)^  , 

where  q^  is  the  radius  of  the  q-ball.  This  type  of  interpolation  has  been  previously  used  and  discussed  in 
(15)-(16).  An  advantage  of  this  model  over  the  original  QBI  model,  i.e.  E{qu)  =  E{u)S{q  -  qo)  (see  (17)), 
is  the  compatibility  with  £(0)  =  1 . 

After  applying  the  mono-exponential  assumption  and  a  few  more  steps  of  calculations  (see  Appendix  C 
for  details),  the  following  ODF  expression  can  be  derived: 

ODFiz)  =  ^  +  ^  r  Vl  ln(-  In  E{{i))d,/>  ■ 

Finally,  rewriting  the  expression  independently  of  the  choice  of  axes,  the  following  analytical  formula 
can  be  shown  to  hold  for  the  ODF: 

ODF{u)  =  —  +  FRT^^l  ln(-  In  £(«))}  [3] 

4;r  16.?r 

where  FRT  is  the  Funk-Radon  transform  (18),  defined  as 

FRT {/(«)}:=  /(vv)j(|vv|  -l)d^w, 
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with  S(»)  the  Dirac  delta  function. 


The  above  ODF  expression  is  dimensionless  and  intrinsically  normalized,  since  the  integrals  of  the  first 
and  second  terms  over  the  sphere  are  respectively  1  and  0.  This  is  in  contrast  to  the  ODF  formulas  used  in 

the  original  QBI,  i.e.,  —  fRt{e{u)},  and  also  in  (14),  where  an  artificial  normalization  factor  Z  is  needed. 

Additional  fundamental  differences  can  he  observed  in  the  approach  presented  here,  compared  to  (14). 
As  we  demonstrated,  integration  of  the  radial  part  of  the  Laplacian  on  the  plane  always  results  in  a 
constant  and  does  not  require  any  model  for  the  diffusion  signal.  Yet,  (14)  uses  the  Bessel  approximation 
of  the  Dirac  delta  function  which  yields  a  variable  (sometimes  even  negative)  term.  As  for  the  integral  of 
the  tangential  term  of  the  Laplacian,  we  use  the  exponential  model  that  is  in  particular  consistent  with 
£:(o)  =  1 ,  in  contrast  to  (14)  that  assumes  the  tangential  term  of  the  Laplacian  to  be  zero  outside  the  q-ball 
(Bessel  approximation  again),  leading  to  an  expression  rather  similar  to  Laplace-Beltrami  post-processing 
sharpening  (19).  A  major  disadvantage  of  approximating  the  Dirac  delta  function  with  a  Bessel  function 
while  considering  the  factor  is  that,  unlike  for  P{f)  which  is  typically  concentrated  near  the  origin, 

the  projection  of  p(r)|fp  may  have  its  highest  values  at  a  certain  positive  radius  coinciding  with  the  side 
lobes  of  the  Bessel  function,  reducing  the  accuracy  of  the  approximation. 

From  Eq.  [3],  it  can  be  seen  that  the  essential  quantity  used  in  computing  the  ODF  from  the  raw  data  is 
l{e)~  ln(-ln.E),  which  is  plotted  along  with  the  absolute  value  of  its  derivative  with  respect  to  £,  in 

Fig.  3.  The  behavior  of  this  quantity  is  almost  linear  for  the  values  of  the  signal  close  to  =0.368, 
which  makes  the  reconstructed  ODFs  corresponding  to  such  signals  similar  to  those  obtained  by  the 
original  QBI  (3)  with  Laplace-Beltrami  sharpening.  Nevertheless,  becomes  greatly  non-linear  as 

the  range  of  the  signal  values  approaches  0  or  1  (for  e.g.  by  acquiring  the  data  at  respectively  higher  and 
lower  b-values),  and  the  advantages  of  Eq.  [3]  become  obvious,  particularly  in  resolving  fiber  crossings 
(see  the  “Results  and  Discussions”  section). 


Implementation 

Our  implementation  of  the  ODE  reconstruction  makes  use  of  the  spherical  harmonic  (SH)  basis,  Tj'"(m), 
which  is  common  for  the  analysis  of  HARDI  data.  The  steps  taken  here  to  numerically  compute  Eq.  [3] 
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are  similar  to  those  described  in  (5).  Particularly,  we  use  the  real  and  symmetric  modified  SH  basis 
introduced  in  (5),  where  SH  functions  are  indexed  by  a  single  parameter  j  =  k{k  +  \)/2  +  m  +  \,  with 
corresponding  kj  and  nij,  as  follows: 


V2Re{7,7l 

V2Im{F“'}, 


-  kj  <mj<0 

Mj  =  0 

0  <mj  <  kj 


We  adopt  a  minimum  least  square  scheme  to  compute  the  modified  SH  coefficients  c  j  of  the  double 
logarithm  of  the  signal,  such  that 

ln(- \nE{u))  ~  ^ CjYj {u),  [4] 

j=i 

where  R-[l  + 1)(/  +  2)/2 ,  with  I  being  the  order  of  the  SH  basis  (we  chose  /  =  4  throughout  our 
experiments).  Next,  since  the  SH  elements  are  eigenfunctions  of  the  Laplace-Beltrami  operator,  we 
compute  ln(-ln£(M))  by  multiplying  the  coefficients  Cj  by  their  corresponding  eigenvalues, 
-kj{kj+l)-  Then,  as  suggested  in  (5),  the  Funk-Radon  transform  is  computed  by  multiplying  the 
coefficients  by  2,^^.  (o),  where  P^(*)  is  the  Legendre  polynomial  of  degree  k,  with 

p  (o)  =  (- 1)2  ^  ^  3  X  ■  •  •  X  (k — ^  ^  Finally,  given  that  v  («)  =  — 1 — ,  the  SH  coefficients  of  the 

2x4x---xk  2-sf^ 

ODF  are  derived  as 


^2 


1 

2^ 

1  .  .h  lx3x---x(k^.  -i-lj 
^  "  2x4x---x{k.  -2) 


7  =  1 

;>i 


By  taking  advantage  of  the  SH  framework,  this  implementation  of  the  proposed  technique  for  the  true 
ODF  is  as  straightforward  as  the  one  introduced  in  (5)  for  the  original  QBI  ODF  formula.  Additionally, 
neither  normalization,  nor  sharpening  is  required  with  this  technique.  This  work  was  recently  extended  in 
(20)  to  impose  positivity  constraint  and  spatial  regularity  if  further  desired. 
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Regularization 


As  mentioned  before,  the  essential  quantity  used  in  computing  the  ODF  from  the  raw  data  is 
l{e)~  ln(-ln£),  plotted  in  Fig.  3.  Hence,  if  there  is  a  relatively  constant  error  A£  in  the  diffusion  data, 
the  error  introduced  in  the  computed  ODF  will  be  proportional  to  the  derivative  of  l{e\. 


al(f)  = 


A£ 

E\nE 


As  Fig.  3  suggests,  this  quantity  becomes  unstable  for  values  of  E  close  to  0  and  1,  which 
subsequently  amplifies  the  error  in  the  diffusion  data.  To  overcome  this  problem,  we  propose  using  a 
flexible  threshold  on  the  diffusion  data  in  order  to  keep  their  values  away  from  the  unstable  regions  of 
[O,  (5j]  and  [l-^2,l],  where  the  thresholds  and  S2  are  manually  defined.  To  perform  this  operation,  we 
use  the  following  function  /(f),  plotted  in  Fig.  4  for  cJj  =  =  0.15 : 


2  ’ 

S,  E^ 

—  + - , 

2  2(^1 


£  <0 
0<E<S^ 
Si<E<l-S2 

1-4  <£<1 
1<£ 


Conversely,  the  ODF  is  most  stable  to  noise  when  Al{e)  is  minimum,  which  is  achieved  for 
E  -  e~^  ~  0.368 .  This  gives  us  a  clue  on  how  to  choose  an  optimum  b-value  in  data  acquisition. 
Particularly,  in  the  mono-exponential  model,  since  E  -  where  the  Apparent  Diffusion  Coefficient 

(ADC)  is  assumed  independent  of  the  b-value,  the  optimum  b-value,  b* ,  is  obtained  as 


{ADC)  ’ 


where  {AD(^  is  the  mean  ADC  in  the  region  of  interest.  Note  that  this  result  holds  only  in  the  simple 
model  which  assumes  both  ADC  and  AE  to  be  independent  of  the  b-value. 
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Extension  to  Multiple  q-Shells 


Multi-Exponential  Model 

We  have  so  far  employed  the  proposed  technique  to  compute  the  ODF  from  a  single  q-shell.  However,  if 
diffusion  data  are  available  on  multiple  q-shells,  this  technique  can  he  applied  to  reconstruct  the  ODF 
while  exploiting  the  information  from  all  the  q-shells.  With  more  available  data,  richer  models  become 
practical  and  appealing.  Here  we  consider  the  following  radial  multi-exponential  model  (see  (16), (21)), 

"  2 

E[qu)  =  Y,K{u)oCkW  ’ 

k=\ 


with  the  constraints 


0  <  <  1, 


z 


[5] 


where  Eq.  [5]  comes  from  the  fact  that  £(0)  =  1 .  Once  the  values  of  and  are  estimated  (see  the 

“Parameter  Estimation”  subsection),  they  can  be  used  in  the  following  more  general  ODE  expression, 
which  is  derived  in  details  in  Appendix  D: 


ODF  {u)  =  ^  +  ^  FFTlvlf^X,  (M)ln(-  In  a,  (m)) 
Att  16;r  rrf 


The  implementation  is  quite  similar  to  what  we  explained  in  the  “Implementation”  section,  with  Eq.  [4] 
being  the  only  difference,  as  it  now  writes 


Z  (M)ln(-  In  a,  (m))  -  c/.(m)  ■ 

k=l  j=l 

In  addition,  the  function  /(e)  introduced  in  the  “Regularization”  section  can  be  applied  to  s,  to 
reduce  the  effect  of  the  noise. 
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Parameter  Estimation 


In  order  to  approximate  the  diffusion  signal  in  a  direction  m  by  a  weighted  sum  of  N  exponentials,  we 
need  to  estimate  the  2N  parameters  \{u)  and  a^.{u)  ,  for  k  =  l,...,N .  To  estimate  the  parameters,  at 
least  2N  - 1  independent  equations  -  besides  Eq.  [5]  -  are  required,  which  can  be  obtained  from  the 
HARDI  signals  measured  on  M  q-balls,  for  M  >  2N  - 1 ,  as  follows: 

^eXu) 

k=l 

where  eXu)'=  E{qfi)  and  q.  corresponds  to  the  i'^  q-ball.  Parameterizing  the  problem  in  terms  of  b- 
values,  bj  =  mX ,  and  choosing  the  physical  units  such  that  the  diffusion  time  becomes  r  =  1 ,  we  obtain 

k=l 

/  =  1,...,M 

Numerical  optimization  approaches  such  as  the  trust  region  algorithm,  (22),  may  be  employed  to  solve 
this  non-linear  system  in  the  most  general  case.  Here,  however,  we  discuss  two  special  cases  (one  familiar 
and  one  new)  with  analytical  solutions.  We  continue  this  subsection  considering  a  fixed  direction,  and 
therefore  drop  the  notation  («). 

The  mono-exponential  assumption  ( =  1)  requires  measurement  on  at  least  M  =  1  q-ball.  M  =  1  leads 
to  Aj  =  1  and  ■  As  it  is  shown  in  Appendix  D,  a/  can  also  be  a  solution  with  any  constant  y. 

Therefore,  choosing  y-bi  results  in  the  solution  ai=  p,  which  is  consistent  with  what  we  already 
derived  (Eq.  [3]).  Eurthermore,  if  measured  values  are  provided  on  more  than  one  q-shells  and  the  mono¬ 
exponential  model  is  still  desired,  then  the  assumption  in  this  model  (ADC  being  independent  of  the  b- 
value)  suggests  that  the  best  exponential  can  be  fitted  by  computing  the  average  ADC  across  all  the  q- 
balls. 

Another  practical  case  of  great  interest  arises  when  we  consider  the  aforementioned  richer  bi¬ 
exponential  model  (N  =  2,  see  for  example  (23)-(24))  to  reconstruct  the  DDEs  from  (at  least)  M  -  3  q- 
shells.  Eor  M  -  3  ,  the  following  system  of  equations  holds  for  each  direction: 
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=  £2 

/la''’+(l-/lK’  =£3 

0<a,^,A<l 

In  general,  this  set  of  equations  can  be  solved  numerically.  Nevertheless,  an  analytical  solution  can  be 
derived  for  the  particular  and  reasonable  case  when  the  sequence  Q,b^,b2,b^  is  an  arithmetic  progress  (the 

sequence  x,  is  an  arithmetic  progress  if  is  constant).  We  describe  this  solution  here,  along  with 

some  regularization  that  guarantees  the  parameters  to  remain  within  the  correct  range. 

Without  loss  of  generality,  let  us  assume  a>  {3 ,  and  also  choose  the  physical  units  such  that  b^-\, 
b^ -2,  and  =  3 .  Then, 


Aa  +  {l-A)jB  =  E^ 
Aa""  +  {l-  ^  E2 

=£3 

0<a,^,A<l 


We  first  define  and  calculate  the  following  two  quantities: 


.4:=^ 

2 

2 


£3  -  EE- 

2fe-k') 

If  e-^AT 


F  F  -F  ^ 
E  -E^ 

1^2 


The  parameters  are  afterward  computed  as  follows: 

a  -  A  +  B 
E^A-B 

A^U^ 

2  2B 

However,  we  still  need  to  ensure  that  they  are  real  and  in  the  correct  ranges.  One  can  verify  that  these 
conditions  are  satisfied  by  enforcing  the  following  constraints: 
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0  <  <l 

<  E, 

E2  <  E^E^ 

£3  -  E^E^  <  £2  -  E^  +  £i£3  -  E2 

Thus,  we  can  obtain  the  optimal  values  of  a ,  (5,  and  X,  by  initially  projecting  £.  s  onto  the  subspace 

defined  by  the  above  inequalities,  and  then  computing  the  parameters.  Note  that  such  projection  is  usually 
necessary,  because  the  bi-exponential  model  may  not  be  fully  accurate  and  the  data  may  be  noisy. 
Furthermore,  using  a  small  separating  margin  of  =  0.01  ~  0.1  in  the  inequalities  makes  the  ODFs  more 
stable  in  practice. 


RESULTS  AND  DISCUSSIONS 

Results  from  Single  q-Shell 

To  validate  our  approach,  we  first  show  results  using  artificial  data.  We  simulated  fiber  crossing  by 
generating  diffusion  images  from  the  sum  of  two  exponentials,  £(m)  =  where  is  a 

diagonal  matrix  with  diagonal  entries  (9,  2,  2),  and  is  rotated  about  the  y-axis  by  a  varying  angle, 

producing  diffusion  values  smaller  than  0.15  (corresponding  to  rather  high  b-values).  The  ODFs  were 
reconstructed  in  the  fourth  order  SH  basis  from  76  diffusion  directions,  uniformly  sampled  on  the  sphere. 
The  results  are  shown  in  Fig.  5,  for  three  different  methods:  our  proposed  framework,  the  original 
(standard)  QBI,  and  the  original  QBI  followed  by  Laplace-Beltrami  sharpening,  (l-TV^)  (see  (19)),  with 

parameter  2=0.15  (chosen  to  produce  the  optimal  results).  As  can  be  seen,  our  method  resolves  the 
crossings  starting  at  about  45°,  compared  to  about  60°  by  the  other  two  methods.  We  also  verified  this 
using  the  dip  test  (25)  -  a  measure  of  multimodality  in  a  distribution  -  on  the  reconstructed  ODFs  from 
the  same  synthetic  diffusion  signals,  with  Rician  noise.  As  can  be  observed  in  Fig.  6,  the  two  modes  of 
the  ODFs  with  high  SNR  are  distinguished  at  a  crossing  angle  which  is  about  15°  smaller  with  the 
proposed  reconstruction  method,  compared  to  the  two  other  techniques.  As  expected,  this  difference 
becomes  less  marked  as  the  noise  increases. 

We  also  tested  our  method  on  three  real  HARDI  datasets;  first  on  the  physical  phantom  in  (26),  which 
was  constructed  from  excised  rat  spinal  cords  and  designed  to  have  crossing  tracts  (90  diffusion  images  at 
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b=1300  s/mm^),  and  then  on  human  hrain  data  (27)  (200  diffusion  images  at  h=3000  s/mm^).  (For  the 
third  real  dataset,  set  the  “Results  from  Multiple  q-Shells”  section.)  The  ODFs  were  reconstructed  with 
the  fourth  order  SH  basis  using  three  approaches:  our  proposed  method,  the  original  (standard)  QBI,  and 
the  original  QBI  followed  by  Laplace-Beltrami  sharpening  with  parameters  0.5  for  the  rat  data  and  0.8  for 
the  brain  data.  Results  are  superimposed  on  the  generalized  fractional  anisotropy  (GFA)  map  and 
presented  in  Fig.  7.  (Note  that  the  ODFs  are  shown  as  they  are;  no  min-max  normalization  is  used  in  any 
of  the  figures.)  Our  method  (left)  produces  sharper  and  more  accurate  ODFs  than  the  original  QBI 
(middle).  In  addition,  although  sharpening  (right)  enhances  the  original  QBI  ODFs  considerably  in 
anisotropic  tissue,  it  causes  significant  instability  in  isotropic  regions  (e.g.  the  background  of  the  rat 
phantom  and  the  CSF  in  the  human  brain  data),  in  contrast  to  our  technique  which  preserves  isotropy 
fairly  well.  For  the  human  brain  dataset,  we  focus  on  the  region  of  the  centrum  semiovale,  where  three 
major  fiber  bundles  intersect:  the  internal  capsule  (IC)/corona  radiata  (CR),  the  radiations  of  the  corpus 
callosum  (CC),  and  the  superior  longitudinal  fasciculus  (SLF). 


Results  from  Multiple  q-Shells 

To  demonstrate  the  advantages  of  exploiting  multiple  q-shells  in  QBI,  we  first  show  experimental  results 
on  a  synthetic  example  which  consists  of  large  diffusion  values  in  two  orthogonal  directions.  We 

synthesized  diffusion  images  by  sampling  the  sum  of  two  exponentials,  =  ||sin^|''  +  |cos(Zi|‘^ 

on  seven  q-shells  {b  =  =l,2,...,l)  and  in  76  directions,  uniformly  distributed  on  the  sphere.  Figure  8 

illustrates  the  ODFs  reconstructed  from  single  q-shells  for  different  b-values,  and  from  three  q-shells  with 
both  mono-exponential  and  bi-exponential  models.  As  can  be  observed,  for  the  data  acquired  at  low  b- 
values  (fo  =  1,2,3),  the  proposed  bi-exponential  model  using  three  q-shells  is  the  only  method  correctly 
resolving  the  horizontal  and  vertical  ODF  peaks,  corresponding  to  the  strong  ADC  values  in  those 
directions  (<^  =  0°,  90°,  180°,  270°).  It  should  be  noted,  however,  that  the  drawback  of  such  a  more  general 
model  is  its  lesser  robustness  to  noise,  as  low  order  models  are  often  more  robust  (e.g.,  computing  the 
average  of  a  signal  is  more  robust  than  estimating  the  actual  signal).  Dark  red  represents  negative  values. 
These  values  do  not  appear  often  in  general,  nonetheless,  a  possible  formal  approach  to  handle  them  can 
be  found  at  (20). 

We  also  tested  our  method  on  the  real  HARDI  dataset  initially  introduced  in  (28).  An  anesthetized 
Macaca  mulatta  monkey  was  scanned  using  a  7T  MR  scanner  (Siemens)  equipped  with  a  head  gradient 
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coil  (80mT/m  G-maximum,  200mT/m/ms)  with  a  diffusion  weighted  spin-echo  EPI  sequence.  Diffusion 
images  were  acquired  (twice  during  the  same  session,  and  then  averaged)  over  100  directions  uniformly 
distributed  on  the  sphere.  We  used  three  h-values  of  1000,  2000,  and  3000  s/mm^,  TR/TE  of  4600/65  ms, 
and  a  voxel  size  of  1x1x1  mm^.  The  proposed  method  was  used  to  reconstruct  the  ODEs  from  the  three  q- 
shells  using  both  bi-exponential  and  mono-exponential  methods,  and  also  from  the  single  q-shells 
individually.  Eigure  9  depicts  the  results  on  a  coronal  slice  through  the  centrum  semiovale  area, 
superimposed  on  the  fractional  anisotropy  (EA)  map.  (Eor  comparison,  one  of  the  sub-figures  shows 
results  by  the  original  QBE)  Note  how  using  the  bi-exponential  method  allows  for  more  clear  recovery  of 
certain  fiber  bundles,  such  as  callosal  radiations  and  corticospinal  tract,  and  better  resolution  of  crossing 
areas  (see  outlined  regions  in  Eig.  9). 
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APPENDIX  A 


Fourier  Transform  of  P(r)|rp 

From  the  Fourier  analysis,  we  know  that  if  E{q)  is  the  Fourier  transform  function  of  P{f),  then: 

F{xP{f)]^i^E{q) 

E{x^P{fp-^E{q) 

dq^ 

where  f{»}  is  the  Fourier  transform  operator.  By  writing  the  second  equation  for  y  and  z  and  summing 
them,  we  will  get: 

E{r^P(r)}=-V^E(q) 

APPENDIX  B 

Computing  the  Radial  Projection  of  a  Symmetric  Function  in  the  Fourier  Domain 

Let  /■:]R^  — he  a  symmetric  function  with  the  3D  Fourier  transform  function  f(q),  and  u  he  an 
arbitrary  unit  vector.  We  will  show  that  [°°f(ru)dr=  ^  ff  f{q)d^q,  where  u'‘~  is  the  plane 
perpendicular  to  u . 

Without  loss  of  generality,  we  choose  our  coordinates  such  that  z-u,  thus  making  the  qx-qy  plane. 
We  first  rewrite  the  expression  as  a  volume  integral  over  the  entire  space,  with  the  help  of  Dirac  delta 
functions: 

Jo  =  £  f{0,0,z)dz  =  ^j\j^j{x,y,z)S{x)S{y)dxdydz’ 
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where  the  factor  —  is  required  because  we  need  the  integral  only  on  the  positive  half  of  the  z-axis,  and  the 
2 

function  is  symmetric.  Let  us  define  g{x,y,z)'=  S{x)S{y)-  For  the  two  functions  with 

Fourier  transform  functions  f{q)  and  g{q),  Parseval’s  theorem  states  that 


z)  dxdydz  =  {q^,q^,q^)  f  {q^,q^,q^)  dqjq^dq^  ■ 

(zTr) 


Computing  g[q^^q^^qj=  2;rS{qJ  and  replacing  it  in  the  above  equations,  leads  to 


£  /  (^2)  dr  =  3'>  z)  g  {x,  y,z)dxdydz 

1 

2{27rf 

"^U-^Hd.’dy,0)dqJq^ 

As  can  be  seen,  the  integral  is  taken  on  the  qx-qy  plane,  which  is  u'‘~ .  This  completes  the  proof. 


E  3  ( q, )  dq^dq/q^ 


APPENDIX  C 

Incorporating  the  Mono-Exponential  Model  in  the  ODF  Formula 

We  will  show  here  that  by  assuming  the  mono-exponential  model,  E{q)  =  e{u)^  ,  we  have: 

[-VlE{q)dqd^  =  \n(-lnE{u))d^, 

q  L 

TT 

while  ^  =  —  is  kept  constant  in  the  integration.  We  begin  by  proving  a  lemma: 

2 

Lemma:  For  a  continuous  and  differentiable  function  / {6,(p) :  5”^  — >  R  with  5^  being  the  unit  sphere, 
we  have: 


where  the  subscript  indicates  the  partial  derivative. 

Proof:  We  use  the  following  expansion  for  Laplace-Beltrami  operator: 
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v,V  =  - 


1  d 


sm 


0d0 


df),  1  aV 


sin6>—  +-  ,  ^  , 

30  J  sin^0d^ 

1 


^cot0.f,  +  f,,+—^f 
sm  0 


^bf  i  ^  =  fee  +  L, 


\0_^  J  66  '  J 
~2 


The  integral  of  the  second  term  is  zero,  because  of  the  periodicity  of  : 

Therefore  the  only  remaining  term  in  the  integral  is  f^g ,  which  completes  the  proof  of  the  lemma. 


We  now  change  the  orders  of  the  integrals  twice,  while  using  the  lemma  in  between: 

f  r  = r  vr^*®*^*^** 

q  q 

_  r  ^  3  p  ^ 

Jo  q  36»Jo  ^ 


[6] 


Next,  we  compute  the  radial  integral: 


We  know  that  0  <  e{u)  <  1  ^  In  e{u)  <  0 ,  so  the  above  expression  vanishes  as  q  ^  oo .  Also,  since  for  a 

negative  function  f{0)  we  have  =  [ln|/((9)|]  =  [ln(-/(^))]g>  the  above  integral  simplifies  as: 

/w)  ^ 
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[  -Ee{q)dq  =  -^[ln(-ln£(M))]^ 
q  2 

Substituting  in  Eq.  [6]: 

=  -|£"[ln(-ln£(M))L# 

=  -irV^ln{-ln£(«)W 

We  completed  the  proof  by  reusing  the  lemma  in  the  last  step. 


APPENDIX  D 


Incorporating  the  Multi-Exponential  Model  in  the  ODF  Formula 


By  assuming  the  multi-exponential  model,  E{qu)  =  '^\{u)ai^{uY  >  we  will  show  that: 


C  r-ylE{q)dqd(t>  =  (M)ln(-  In  a,  {u))d^ , 

q  ^  k=i 


TT 

while  ^  =  —  is  kept  constant  in  the  integration.  The  ODF  will  then  be  derived  by  replacing  the  above 
2 

expression  in  Eq.  [2]. 

The  proof  is  an  extension  of  Appendix  C.  We  first  introduce  the  new  non-negative  variable 
5^(m) :=  -lna^(M)  which  yields  E{qu)  =  ■  For  simplicity,  here  we  drop  the  notation  («).  We 

k=l 

then  continue  from  Eq.  [6]  and  compute  the  radial  integral 


q\k=i 

N 


dq 


[7] 


k=l  q  k=l 


dq 


The  first  integral  is  computed  the  same  way  as  in  Appendix  C: 
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k=l  k=l 


N  . 

f"i  2s,  10 

^  i=i 


Regarding  the  second  integral  of  Eq.  [7],  let  us 


define  l{s):^r-f^\^e-^^^^dq 

^  k=l 


with  s  the  vector  of 


s,  s,  and  derive  it  with  respect  to  s,  '■ 


Then,  we  can  see  that: 


ds. 


'ke'^ 


A  2 


dq 


2s, 

A. 


■ke 

2s 


ds. 


k 

^  A  ^ 

^Ins 
2 


^  k=l 


where  C  is  independent  of  s  .  By  evaluating  the  function  for  s  =1  (vector  of  all  Is),  we  obtain  C  =  /(l), 
which  we  then  compute  using  the  original  definition  of  l{s): 


l{s) 


^  k=\ 
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Therefore, 


We  now  insert  the  values  of  the  two  integrals  in  Eq.  [7]: 

j  —Egdq--— 

q  ^  k=i  ^  t=i 

\  (  N  A 

=  “9  ZAln^ 

^Vr=i  Je 

I  7  iv  ^ 

=  ZAln(-ln«J 

^Vr=i  Je 

Finally,  substituting  in  Eq.  [6] : 

^-'flE{rq]ded^ 


We  completed  the  proof  hy  using  the  lemma  introduced  in  Appendix  C. 


An  interesting  observation  is  that  if  ^^(m)  is  a  set  of  estimated  parameters  in  the  multi-exponential 
model,  then  for  any  constant  y,  the  set  (X^{uY  results  in  the  same  ODE: 

=  {uYn{~ylna,  («)) 

k=l  k=l 

^^lj^^{u)[in{-lna,{u))  +  lny~j 

k=l 

^Vlf^A,{u)ln{-lna,{u))  +  Vl(lnrf^A,{u) 

k=i  V  k=i 

-ylt,A,{{i)\n{-\na,{{l))  +  Vllny 

k=l 

k=l 

This  is  expected,  since  the  ODF  is  dimensionless  and  should  not  depend  on  the  physical  units  of  q. 


d<f> 
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(a)  (b) 

Fig.  1.  Radial  integration  of  the  PDF,  (a)  in  a  cone  of  constant  solid  angle  (i.e.,  the  factor  is  considered),  and  (b)  by  linear 
projection  (i.e.,  inaccurately  without  the  factor  r^as  commonly  done  in  the  HARDI  literature  before  this  work). 
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Fig.  2.  DTI  example  of  ODF  reconstruction  (with  (10,  5,  Ij  as  the  diagonal  entries  of  the  tensor),  shown  from  two  view  angles,  (left) 
with  the  factor  r^,  (right)  without  the  factor  and  after  normalization.  Note  how  less  sharp  the  latter  is  and  how  incompletely  it 
represents  the  true  structure  of  the  ODF. 
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Fig.  3.  Behavior  o/ln(—  InF)  (left)  and  the  absolute  value  of  its  derivative  (right)  with  respect  to  E.  Note  how  unstable  they  are  for  E 
close  to  0  or  1. 
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Fig.  4.  The  regularization  function  used  for  the  diffusion  signal  to  avoid  the  unstable  regions  (blue  curve). 
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Fig.  5.  Experimental  results  on  synthetic  data  with  fiber  crossing,  using:  (top)  our  proposed  technique,  (middle)  original  QBI  after 
normalization,  and  (bottom)  original  QBI  with  Laplace-Beltrami  sharpening. 
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Fig.  6.  Results  of  the  dip  test  on  the  same  ODFs  as  in  Fig.  5,  with  various  noise  levels. 
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Fig.  7.  Reconstructed  ODFs  from  rat  spinal  cord  phantom  (top)  and  human  brain  (bottom),  shown  on  the  GFA  map,  using:  (left)  our 
proposed  technique,  ( middle )  original  QBI  after  normalization,  and  ( right)  original  QBI  with  Laplace-Beltrami  sharpening. 
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mono-exp.  bi-exp. 

b  =  l  b  =  2  b  =  3  b  =  4  b  =  5  b  =  6  b  =  7  b=  1,2,3  b  =  1,2,3 

Fig.  8.  Results  of  the  ODF  reconstruction  on  synthetic  data.  Note  how  the  bi-exponential  model 
correctly  resolves  the  maxima  of  the  ODF  from  low  b-values. 
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Fig.  9.  Reconstructed  ODFs  from  the  real  brain  data,  shown  on  the  FA  map.  The  bi-exponential 
model  ODFs  (top,  left)  have  been  scaled  down  1.5  times  for  better  comparison.  All  the  ODFs 
except  those  in  (top,  right)  have  been  reconstructed  considering  the  factor  r^.  Note  how  the  bi¬ 
exponential  model  for  diffusion  improves  the  resolution  of  fiber  crossings,  compared  to  the 
mono-exponential  (constant  ADC)  model. 
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